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A  MULTI-DOMAIN  SPECTRAL  METHOD  FOR  SUPERSONIC  REACTIVE  FLOWS 


WAI-SUN  DON,  DAVID  GOTTLIEB  k  JAE-IIUN  JUNG 


Abstract.  This  paper  has  a  dual  purpose:  it  presents  a  multidomain  Chebyshev  method  for  the  solu¬ 
tion  of  the  two-dimensional  reactive  compressible  Navier-Stokes  equations,  and  it  reports  the  results  of  the 
application  of  this  code  to  the  numerical  simulations  of  high  Mach  number  reactive  flows  in  recessed  cavity. 
The  computational  method  utilizes  newly  derived  interface  boundary  conditions  as  well  as  an  adaptive  fil¬ 
tering  technique  to  stabilize  the  computations.  The  results  of  the  simulations  are  relevant  to  recessed  cavity 
flameholders. 

Key  words,  multi-domain  spectral  method,  penalty  interface  conditions,  supersonic  combustor,  recessed 
cavity  flame-holder,  compressible  Navier-Stokes  equations 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  The  efficacy  of  spectral  methods  for  the  numerical  solution  of  highly  supersonic, 
reactive  flows  had  been  previously  reported  in  the  literature.  Don  and  Gottlieb  [7,  8]  simulated  interactions 
of  shock  waves  with  hydrogen  jets  and  obtained  results  showing  the  rich  dynamics  of  the  mixing  process  as 
well  as  the  very  complex  shock  structures.  Don  and  Quillen  [9]  studied  the  interaction  of  a  planar  shock 
with  a  cylindrical  volume  of  a  light  gas  and  showed  that  the  spectral  methods  used  gave  good  results  for  the 
flows  with  the  shocks  and  complicated  non-linear  behaviors.  In  fact  the  results  compared  favorably  to  ENO 
schemes. 

The  methods  reported  above  were  based  on  Chebyshev  techniques  in  one  domain.  In  order  to  extend 
the  utility  of  spectral  methods  to  complex  domains,  multidomain  techniques  have  to  be  considered.  The 
main  issue  here  is  the  stable  imposition  of  the  interface  boundary  conditions,  and  in  this  paper  we  consider 
mainly  the  penalty  method,  introduced  for  hyperbolic  equations  by  Funaro  and  Gottlieb  [10,  11]. 

There  is  an  extensive  literature  on  the  subject:  Hesthaven  [13,  14,  15]  applied  penalty  BC  for  Chebyshev 
multidomain  methods  using  the  characteristic  variables.  Carpenter  et.  ah  [4,  17,  18]  used  it  in  conjunction 
with  compact  finite  difference  schemes,  going  from  a  scalar  model  equation  to  the  full  N-S  equations  in 
general  coordinate  systems.  Carpenter,  Gottlieb  and  Shu  [5]  demonstrated  the  conservation  properties  of 
the  Legendre  multidomain  techniques. 

In  the  current  work  we  follow  the  same  methodology  but  in  the  context  of  supersonic  combustion. 
We  formulate  the  stable  interface  conditions  based  on  the  penalty  method  in  a  conservative  form  for  both 
Euler  and  Navier-Stokes  equations  in  two  dimensional  Cartesian  coordinates.  We  derive  stability  conditions, 
independent  on  the  local  flow  properties,  for  the  penalty  parameters  for  the  Legendre  spectral  method. 
We  also  present  here  a  new  adaptive  filtering  technique  that  stabilize  the  spectral  scheme  when  applied  to 
supersonic  reactive  flows. 

Implementing  this  method,  we  consider  supersonic  combustion  problems  in  recessed  cavities  in  order  to 
establish  the  efficacy  of  recessed  cavity  flame-holders. 

*  Brown  University,  Division  of  Applied  Mathematics,  182  George  Street,  Providence,  RI  02912  (E-mail:  wsdon,  dig, 
jung@cfm.brown.edu)  This  work  was  performed  under  AFOSR  grant  no.  F49620-02-1-0113  and  DOE  grant  no.  DE-FG02- 
96ER25346.  This  research  was  supported  by  the  National  Aeronautics  and  Space  Administration  while  the  second  author  was 
in  residence  at  ICASE,  NASA  Langley  Research  Center,  Hampton,  VA  23681. 
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We  consider  two  different  cases;  (1)  Non-reactive  flows  with  two  chemical  species  and  (2 )  Reactive  flows 
with  four  chemical  species. 

Recessed  cavities  provide  a  high  temperature,  low  speed  recirculating  region  that  can  support  the  pro¬ 
duction  of  radicals  created  during  chemical  reactions.  This  stable  and  efficient  flame-holding  performance 
by  the  cavity  is  achieved  by  generating  a  recirculation  region  inside  the  cavity  where  a  hot  pool  of  radicals 
forms  resulting  in  reducing  the  induction  time  and  thus  obtaining  the  auto-ignition  [2,  22].  Experiments  have 
shown  that  such  efficiency  depends  on  the  geometry  of  the  cavity  such  as  the  degree  of  the  slantness  of  the 
aft  wall  and  the  length  to  depth  ratio  of  cavity  L/D.  Thus  one  can  optimize  the  flame-holding  performance 
by  properly  adjusting  the  geometrical  parameters  of  the  cavity  flame-holder  system  for  a  given  supersonic 
flight  regime.  There  are  two  major  issues  of  such  cavity  flame-holder  system  that  need  to  be  investigated  ; 
(1)  What  is  the  optimal  angle  of  the  aft  wall  for  a  given  L/D?  and  (2) How  does  the  fuel  injection  interact 
with  cavity  flows ?  An  answer  to  these  questions  require  both  a  comprehensive  laboratory  and  numerical 
experiments. 

There  have  been  previous  numerical  studies  on  these  questions,  many  of  them  rely  on  the  turbulence 
models.  Rizzetta  [19]  used  a  modification  of  the  Baldwin  Lomax  algebraic  turbulence  model.  Davis  and 
Bowersox  [6]  also  used  Baldwin-Lomax  model.  Zhang  et.al.  [23]  used  Wilcox  k  —  uj  turbulence  model.  Baurle 
and  Gruber  [3]  used  the  Menter  model.  Although  the  use  of  the  turbulence  models  can  make  it  possible 
to  handle  the  compressible  supersonic  shear  flows,  the  results  are  quite  model-dependent  as  they  require 
parametric  assumptions.  In  this  work,  we  solve  the  full  compressible  Navier-Stokes  equations  with  chemical 
reactions  without  any  turbulence  model,  using  a  multi-domain  spectral  method. 

Results  of  several  numerical  studies  including  the  present  study  have  shown  that  the  stability  of  the 
recirculation  inside  cavity  is  enhanced  for  the  lower  angle  of  cavity  compared  to  the  rectangular  cavity. 
The  present  study,  however,  gives  more  accurate  and  finer  details  of  the  fields  than  those  done  by  lower 
order  numerical  experiments.  We  show  that  a  stationary  recirculation  region  is  not  formed  inside  the  cavity 
contrary  to  what  the  lower  order  schemes  predict.  A  quantitative  analysis  made  in  this  study  shows  that 
the  lower  angled  wall  of  the  cavity  reduces  the  pressure  fluctuations  significantly  inside  the  cavity  for  the 
non-reactive  flows.  We  obtained  a  similar  result  for  the  reactive  flows  with  the  ignition  of  the  fuel  supplied 
initially  in  the  cavity. 

The  rest  of  this  paper  is  organized  as  follows.  In  section  2  the  governing  equations  are  given.  In  section 
3  we  describe  the  numerical  method  used  in  this  work.  In  this  section  we  present  the  adaptive-filtering 
used  to  remove  the  high  frequency  mode  that  causes  the  instability  due  to  the  non-smoothness  of  the  flow, 
and  we  derive  stable  penalty  interface  conditions.  In  section  4  the  system  of  the  supersonic  recessed  cavity 
combustor  is  described.  In  section  5  the  main  results  of  this  work  are  given  and  discussed. 

2.  The  Governing  Equations.  In  this  work,  we  consider  the  compressible  Navier-Stokes  equations 
in  the  presence  of  the  chemical  reactions.  Since  Hydrogen  is  used  as  a  fuel  in  our  numerical  experiments, 
four  chemical  species  are  considered,  i.e.  H2, 02,  H20  and  N2  with  the  chemical  reaction  between  Hydrogen 
and  Oxygen  gases: 


(2.1) 


2H2  ~f"  O2  2H2O. 


The  two-dimensional  compressible  Navier-Stokes  equations  in  conservative  form  can  be  written  as: 


(2.2) 


dq  dlF  dG_dF±  dG^ 
dt  +  dx  +  dy  dx  +  dy  + 


2 


The  state  vector,  q  and  the  in  viscid  fluxes,  F  and  G  are  given  by 
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Here  p  is  the  density,  u  and  v  are  the  mean  mixture  velocity  components  of  flow,  E  is  the  total  internal 
energy  and  P  is  the  pressure.  The  mass  fraction  vector,  is  f  =  (/i,/2,/3,/4)r  and  the  column  vectors  fu 
and  fv  are  composed  of  the  specific  momentum  of  ith  species 


(2.4)  fU£  —  fi(lL  4*  fvi  —  fi(v  4*  fii)* 

The  velocity  field  (ui,Vi)  of  the  ith  species  is  the  drift  velocity  relative  to  the  mean  mixture  velocity  (a,v) 
and  is  determined  by 

(2-5)  (M)=4-V/i. 

P$c. 

Here  p  is  the  mixture  dynamic  viscosity  to  be  determined  in  (2.11),  and  Sc  is  the  Schmidt  number  which  is 
taken  to  be  0.22.  The  viscous  fluxes,  Fv  and  Gv  are  given  by 


Fv  = 


lyx 
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where  0  =  (0,0,0,0)T,  T  is  the  temperature,  Cp  is  the  mixture  specific  heat  at  constant  pressure,  Pr  is 
the  Prandtl  number  (which  is  taken  to  be  0.72)  for  the  normal  air  and  hi  is  the  specific  enthalpy  of  the  ith 
species  and  given  by 


hj  —  h ®  4* 


ds . 


where  h J  is  the  reference  enthalpy  of  the  ith  species  and  the  specific  heat  of  the  ith  species  at  constant 
pressure,  CPi  is  represented  as  a  fourth-order  polynomial  of  T  (see  [16]).  The  elements  of  the  viscous  stress 
tensor  are  given  by 


(2.7) 


+  A  >  V  duk 
+  6iix2^dx  ’ 

k=l  k 


where  5  is  the  Kronecker  delta  symbol,  and  A  is  the  bulk  viscosity  which  is  taken  to  be  - \p  under  the  Stokes 
hypothesis. 
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The  equation  of  state  is  given  by  the  assumption  of  the  perfect  gas  law 

4 

(2.8)  P  =  pBT  =  RT^pfi/Mi, 

i- 1 

where  R  is  a  mixture  gas  constant  with  the  universal  gas  constant  R  and  M*  is  the  molecular  weights  of  ith 
species.  The  energy  E  is  given  by 

(2.9)  E  =  p  [T  Cp(s)ds  -P+  \p{u- 2  +  v2)  +  Y,  PWl 

Z  i=  1 

where  the  mixture  specific  heat  at  constant  pressure  is  given  by 

4 

(2.10)  Cp  =  Y,CPifi!Mi. 

i- 1 

2.1.  The  chemical  models.  We  use  the  same  models  as  in  [7].  Each  chemical  species  has  different 
dynamical  viscosity  /i*  based  on  Sutherland’s  law  and  we  obtain  the  mixture  viscosity  /i  by  Wilke’s  law  [21], 


=  (Ty~  (Tn  +  SA 

Pot  Vo i)  \  T  +  Si  )  ' 

(2.H)  mfi/Mi 

(l  +  [(M/PjXfi/mHMi/Mj)1*)2 

**  ~  [8(1  +  (Mi/Mj))]*  ‘ 

Here  /ioj,  To,-  and  Si  are  constants.  A  modified  Arrhenius  Law  gives  the  equilibrium  reaction  rate  ke,  the 
forward  reaction  rate  kf  and  the  backward  reaction  rate  kb  as 

ke  =  AeT  exp(4.60517(iJe/r  -  2.915)) 


kf  =  Af  exp  (-Ef/{RT)) 
kb  =  kf/ke, 


where  the  activation  energy  Ee  —  12925,  Ef  =  7200  and  the  frequency  factor  Ae  =  83.006156,  Af  =  5.541  x 

1014. 

The  species  are  ordered  as  follows  :  (H2,02,H20,N2),  and  the  law'  of  mass  action  is  used  to  find  the 
net  rate  of  change  in  concentration  of  ith  species  (7,  by  the  single  reaction  (2.1),  i.e. 


Cx  =  2(A/[H2]2[02]  -  4[H20]2) 
C2  = -(^[H2]2[02]  -  MH2O]2) 
C3  =  2(fc/[H2]2[02]  -  MH2O]2) 


where  [•]  denoted  the  net  rate  of  change  in  concentration. 

Finally,  the  chemical  source  term  C  in  (2.2)  is  given  by 

(2.12)  C  =  (0,0,0, 0,CiM1,C2M2,C3M3, , 

where  C%  is  the  net  rate  of  change  in  concentration  of  ith  species  by  the  reaction. 

In  Appendix  C,  a  table  of  all  the  necessary  coefficients  and  constants  used  for  the  reactive  Navier-Stokes 
equations  with  species  (H2,  CbjEbO,  N2)  are  given. 
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3.  The  Multi-domain  Spectral  Method.  In  this  section  we  describe  the  two  crucial  components  of 
the  Chebyshev  multi  domain  code  used  in  our  work,  i.e.  the  adaptive  filtering  and  the  penalty  method  for 
the  stable  interface  conditions. 

3.1.  The  adaptive  filtering  technique.  It  is  well  known  that  spectral  methods  may  exhibit  insta¬ 
bilities  when  applied  to  nonlinear  equations.  To  stabilize  the  spectral  scheme  in  an  efficient  way  we  use  here 
filters  to  attenuate  the  high  frequency  modes  of  the  function  qN(x,t)  smoothly  to  zero.  Thus  the  filtered 
version  of  a  polynomial  qN  is  given  by: 

N 

(3.1)  <£(M)=  X)  *(k/N)ak(t)ct>k(x), 

k=—N 

where  a*  is  the  transform  coefficient  and  <j>k  is  the  basis  polynomial  of  order  k  (generally  the  Fourier  and 
Chebyshev  polynomials  for  a  periodic  and  non-periodic  function  respectively). 

Following  Vandeven  [20]  we  define  a  filter  function  cr(u)  of  order  p  >  1  as  a  C°°[-~  1, 1]  function  satisfying 

(32)  <7(0)  =  1,  *(±i)=o, 

[‘}  cr^)  (0)  =  0 ,  cr^(±l)  =  0,  j<p 

where  denotes  the  j-th  derivative. 

It  can  be  shown  that  the  filtered  sum  (3.1)  approximates  the  original  function  very  well  away  from  the 
discontinuities.  A  good  example  of  filter  function  is  the  exponential  filter.  It  is  defined  as 

(3.3)  cf(u)  =  exp  (-ct|a>|7)  , 

where  -1  <  a?  =  k/N  <  1,  a  =  -lne,  e  is  the  machine  zero  and  7  is  the  order  of  the  filter. 

The  exponential  filter  offers  the  flexibility  of  changing  the  order  of  the  filter  simply  by  specifying  a 
different  7.  One  does  not  have  to  write  a  different  filter  for  different  order.  Thus  varying  7  with  N  yields 
exponential  accuracy  according  to  [20].  In  the  present  study  the  sixth  order  global  smoothing  (7  =  6)  is 
used.  If  the  order  of  the  filter  7  is  taken  to  be  too  small,  say  7  <  4,  the  method  becomes  too  dissipative. 

In  the  current  application,  the  interaction  of  the  aft  cavity  wall  and  the  strong  vortex  generated  by  the 
shear  layer  flow  over  the  cavity,  creates  large  pressure  variations  near  the  corner  of  the  aft  cavity  wall.  The 
local  sharp  gradient  can  cause  numerical  instability  and  a  heavier  filter  is  needed  to  prevent  the  development 
of  oscillations  in  this  region.  This  heavy  filtering  can  be  used  globally  and  maintain  the  stability  of  the 
scheme,  however  this  dissipates  out  all  fine  scale  structures,  which  is  highly  undesirable  when  the  resolution 
of  fine  scale  structures  is  essential  for  the  understanding  of  the  recessed  cavity  flameholder  systems. 

Since  this  is  a  local  phenomenon,  it  is  enough  to  apply  a  heavy  filter  only  in  points  in  this  region. 
This  Local  Adaptive  Filtering  keeps  the  scheme  stable,  without  dissipating  fine  scale  features  away  from  this 
region.  The  local  adaptive  filtering  is  carried  out  where  conditions  such  as  ql  <  q  <  qu  are  violated.  Here  q 
can  be  the  mass  fraction  of  each  species  fi  and/or  temperature  T  and  ql  and  qu  denote  the  lower  and  upper 
tolerance  limits  of  q.  In  this  work,  a  filtering  of  the  order  7  =  2  or  7  =  3  is  used  to  reduce  the  magnitude  of 
the  oscillations  at  those  points. 

The  results  of  this  work  indicate  that  the  local  adaptive  filtering  is  applied  only  in  a  few  number  (in  the 
range  of  1  to  7)  of  grid  points  around  the  corner  of  the  aft  wall  once  in  a  while. 

3.2.  Stable  interface  conditions.  In  this  paper  we  use  mainly  the  penalty  type  interface  conditions, 
i.e.  the  boundary  conditions  are  imposed  only  in  a  weak  form  [10, 11].  Successful  penalty  interface  conditions 
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were  constructed  based  on  the  characteristics  for  the  Navier-Stokes  equations  in  [13,  14,  15]  and  for  spectral 
method  and  for  high-order  finite  difference  methods  in  [4, 17, 18],  and  a  conservative  form  of  penalty  interface 
conditions  was  proposed  [5]  for  the  Legendre  spectral  method.  Following  the  same  idea  as  those  works,  we 
consider  two  interface  conditions,  i.e. 

1.  The  averaging  method ,  in  which  the  interface  conditions  are  obtained  by  averaging  the  state  vectors 
of  the  two  adjacent  domains,  and 

2.  The  Penalty  method  in  conservative  form  in  which  the  interface  conditions  are  satisfied  only  in  a 
weak  form,  leaving  the  approximations  not  necessarily  continuous  at  the  interfaces. 

In  the  following  sections  we  will  give  the  penalty  interface  conditions  for  the  Euler  and  Navier-Stokes 
equations  and  also  show  that  the  averaging  method  is  a  subset  of  the  penalty  method. 


3.2.1.  Conservative  penalty  interface  conditions.  Consider  equation  2.2  with  the  inviscid  part 
only,  in  the  ^-direction  in  the  interval  —  2  <  x  <  2,  i.e., 


(3.4) 


fi!  +  *£=0 

a  Si 


For  simplicity,  assume  that  we  have  two  domains  in  this  interval  with  the  interface  at  x  =  0,  qJN(x,t) 
denotes  the  numerical  solution  in  the  left  domain  x  <  0  and  q{f(x,t)  in  the  right  domain  x  >  0.  Note  that 
the  numerical  solution  is  composed  of  two  polynomials  of  different  orders.  The  Legendre  spectral  penalty 
method  is  given  by 


(3.5) 


dt  dx 


=  B(qIN(-2,t))  + 

t\Qn(x)  [/+(<7n(0i  t))  -  ))]  + 

T2QN(x)[f-{qIN(0,t ))  -  /"(«M(0,t))], 


dt  dx 


T3<2m(z)[/+(<Zm(0,*))  -  /+(«t(M))]  + 


nQM(x)[f  (qlM(0,t))-f  (qlN(0,t))] 


where  B  is  a  boundary  operator  at  the  end  points,  i.e.,  x  —  ±2  and  lfN  and  I  If  are  the  Legendre  interpolation 
operators  for  the  left  and  right  domains  respectively.  .  The  positive  and  negative  fluxes  /+  and  /"  are  defined 
by 

(3.6)  /±  =  j  SA±S~1dq, 
with 

dF 

(3.7)  A=^=SAS~\ 

The  Jacobian  matrix  A  is  assumed  to  be  symmetric.  A+  and  A“  are  the  diagonal  matrices  composed  of 
positive  and  negative  eigenvalues  of  A  respectively.  Qn(x)  and  Qm{%)  are  polynomials  of  orders  N  and  M 
respectively  such  that  they  are  zero  at  all  the  collocation  points  except  the  interface  points  x  =  0  (for  example 
Qn(x)  =  ?  0  <  x  <  2  where  Tn(x)  is  the  Chebyshev  polynomial  of  degree  N).  The  penalty 

parameters  ti,T2,T3  and  r 4  are  all  constants.  Since  we  are  interested  only  in  the  interface  conditions,  we 
ignore  the  boundary  operator  B  at  x  —  ±2,  Define  the  discrete  scalar  product  (p,g)jv  = 

Ui  is  the  weight  in  the  Gauss-Lobatto-Legendre  quadrature  formula.  With  the  discrete  product,  the  energy 
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E(t)  is  defined  by  E(t)  =  (qif(x,t),qIN(x,t))N  +  (QM{x,t),qjj(x,t))M-  The  stability  conditions  of  penalty 
parameters  are  given  by  the  following  theorem  [5]: 

Theorem  3.1.  The  energy  is  hounded  by  the  initial  energy  of  the  system  if  the  following  conditions  are 
satisfied  ; 

2ujjnti  <  1,  2ujjnt2  >  1,  2^73  <  -1,  2^41 74  >  -1, 

(3.8)  oj’nti  -  ojtfr3  =  1,  w^r2  -  =  1- 


3.2.2.  The  penalty  method  for  the  Euler  Equations.  The  penalty  method  in  the  case  of  the  2-D 
Euler  equation  is  given  by 

du  +  Us* M  +  ej* gw)  _  nMl'S)[f+M  _  /+(„„_)]  + 

(3.9)  T2AQ(x,y)[f~ (qN)  -  f~(qM -)]> 


where  qM~  is  the  state  vector  of  the  adjacent  domain  at  the  interface  of  degree  M,  71,3(72,4)  denotes  71(72) 
and  73  (7-4)  respectively,  ti  and  r<±  (7*3  and  74)  are  the  penalty  parameters  for  the  right(left)  in  ^-direction 
and  top(bottom)  in  y-direction  respectively.  Q(x,y)  is  a.  polynomial  which  vanishes  at  all  of  interior  points 
of  the  domain  and  is  equal  to  1  at  the  four  interfaces.  Note  that  the  boundary  operator  B  does  not  appear 
in  the  scheme.  Let  A  be  the  linearized  Jacobian  matrix  (around  a  state  vector  q0)  of  two  inviscid  fluxes 


(3.10) 


(dF  dG\ 
\dq'dq) 


where  n  =  ( nX7ny )  is  the  unit  outward  normal  vector.  Since  the  matrix  A  is  symmetric,  there  exists  S  such 
that 


(3.11)  A  —  SAS"1, 

where  A  is  a  diagonal  matrix  composed  of  eigenvalues  of  A.  Then  A  =  +  A~~  and  A±  —  SA±S~1.  A±  is 

defined  as  in  previous  section.  Splitting  A  yields 

(3.12)  f±  =  A±q0,. 


where  f±  is  obtained  from  the  linearized  state. 

Remark  1.  Since  ft  =  (nx,ny)  is  taken  to  be  outward  normal  vector ,  the  stability  condition  (3.8)  is  now 
modified  and  given  as 

2ojtnTi  <  1,  2 oj*nT2  >  1,  2cj^r3  <  1,  2wj^r4  >  1, 

(3.13)  +  CJmT4  =  1,  ujsfT‘2  +  VmTs  =  1. 


The  Jacobian  matrix  A  and  its  eigenvalue  matrix  A  are  given  in  Appendix  A. 

For  illustration,  we  consider  the  propagation  of  a  Gaussian  density  peak  at  the  center  of  rectangular 
physical  domains.  The  physical  domain  is  partitioned  with  16  sub-domains.  The  interface  conditions  be¬ 
tween  the  domains  are  imposed  according  to  the  penalty  Euler  equations  as  discussed  above.  Characteristic 
boundary  conditions  are  imposed  at  the  outer  physical  boundaries.  The  results  presented  in  figure  3.1  indi¬ 
cate  that  the  penalty  formulation  works  well.  From  the  numerical  experiments  of  this  problem,  we  observe 
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Fig.  3.1.  The  propagation  of  a  density  peak  with  the  penalty  Euler  equations  with  16  sub-domains:  The  initial  condition 
(left)  and  the  solution  (right)  at  t  —  0.03604ra$  are  given . 

that  reflections  can  be  created  at  the  interface  across  the  adjacent  domains  depending  on  the  choice  of  the 
penalty  parameters.  Thus  proper  choice  of  the  penalty  parameters  should  take  into  account  reflections  from 
the  interfaces.  To  demonstrate  the  above  formulation  for  the  Euler  Equations,  We  will  return  to  this  issue 
in  a  future  paper  1 . 

3.2.3.  The  penalty  method  for  the  Navier- Stokes  Equations.  When  dealing  with  the  Navier 
Stokes  equation,  we  keep  the  penalty  form  for  the  Euler  fluxes  and  add  a  penalty  term  for  the  viscous 
fluxes.  The  stability  of  this  procedure  stems  from  the  fact  that  the  Jacobian  matrices  for  the  full  reactive 
Navier-Stokes  equation  can  be  symmetrized  by  the  same  similarity  transformation  (see  Appendix  B).  Thus 
we  get  the  system: 

Oqn  01  nF  01  nG  _  0InFv  OInGv 

dt  Ox  dy  Ox  dy 

n.3Q{x,y)[f+(qN )  -  f+(qM-)]  + 

T2AQ(x,y)[f~(qN)  -  f~(qM-)]  + 

T6.$Q(x,  y)[ A„  •  qw  -  A„  •  qjvf_]  + 

(3.14)  r5jQ(x,y){Av  •  dqN  -  Av  •  d<\M_]. 

Here  are  same  as  defined  in  the  previous  section  and  the  Jacobian  matrix  vector  A„  is  given  by 

A  fdFu  dGv  \ 

(3  5)  Ur1*’  dqvUV 

and 

(3.16)  q=(q,q),  Oq=(qx,qy), 

where  again  q_  and  <9q_  denote  the  adjacent  domains  state  vectors  and  their  derivatives.  Note  that  the 
penalty  terms  Au  •  dq  does  not  appear  in  [4,  17,  18].  The  penalty  parameters  r5f7  and  are  defined  in 
the  same  way  as  in  the  previous  section.  To  seek  stable  penalty  parameters  we  split  the  inviscid  and  viscous 
fluxes  and  keep  the  stability  conditions  of  71,2,3,4  for  the  inviscid  flux  as  in  Theorem  3.1.  The  stability 
conditions  of  75,7  and  7*6,8  are  given  in  the  following  Theorem  : 

lJ,  H.  Jung,  PhD  thesis,  Div.  of  Applied  Math.,  Brown  University,  2002 
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Theorem  3.2.  The  penalty  method  for  the  Navier-Stokes  equations  (3.14)  is  stable  if  the  penalty  pa¬ 
rameters  Tj  ,  j  =  1...4  are  as  in  Theorem  3.1  and  the  rest  satisfy : 

UnTq  <  0, 

UJNTq  —  COAfTg  =  0, 

1  4  OJNT 5  ~  OJMT7  =  0, 

(3.17)  ( - 1 - ^  ^Mr7  —  2r7  4  4tdjv?6  H - <  0  . 

\umojn)  um 


Proof 

As  in  the  proof  of  Theorem  3.1,  we  assume  that  we  have  two  domains  and  by  multiplying  the  equations 
by  the  state  vectors,  we  get 

(3.18)  <  [Invisdd\  -f  [Viscous], 

2  at 

where  [Inviscid]  and  [Viscous]  denote  the  terms  from  inviscid  and  viscous  parts  of  the  equation  respectively. 
The  conditions  for  rit2  and  7-3,4  given  in  the  Theorem  3.1  assure  that  the  first  term  [Inviscid]  is  negative. 
The  [Viscous]  part  at  the  interface  is  given  by 

JV 

[Viscous]  =  qT  Auq(  -  ^  q’?  A^Uj  - 

i~  0 
M 

qT.Avq'_  -  Yiq'i-T Ai,q'i-Uj  + 

3=  0 

(3.19)  T5uiN[qT Avq'  -  qT Avq'_]  +  TtU}M[qZAvq'_  -  qT.Avq']  + 

TeuN[qTAvq  -  qTA„q _]  +  TgWM[q~Auq-  -  q*Avq], 


where  q '  denotes  the  derivative  of  q  either  in  x  or  y  direction,  u)  is  the  Legendre  weight,  and  Av  is 


(3.20) 


Since  all  the  eigenvalues  of  Au  are  non-negative,  every  term  inside  the  summations  in  the  above  equation  is 
not  negative,  and  we  would  like  to  keep  the  boundary  terms.  Thus  we  get  the  energy  estimate  such  as 


(3.21) 


“£(*)  <  [qTAuq*  -  uNq,T Avql  -  q^Auql_  -  uMql- Auq'_]  4 

T*>u)N[qT Avq*  -  qrAuqt_ ]  4  T?LJAf[q-Avq^  -  q-Avq']  4 
TeUN[qT Auq  -  qT AuqJ]  4  T&WA#[glA„g_  -  ql A„q]  . 


The  RHS  of  (3.21)  can  be  rewritten  as 

(3.22)  RHS  =  utBu , 


where  u  =  (q.q-yq^qL)  and  B  is  given  by 


(3.23) 


'  2 06Au 

-OqAI  -  (7 SAU 

(1  4-  <To)A„ 

\  Al 


-cjqAu  -  ogAl 
2(78  Av 
-<J7  Al 

(“1  4  #7)  Af) 


(1  4  o*f)Av 
-07  Av 

—  2  0JjvAu 

0 


—  (JVyAp  ^ 

(— 1  4  (77 )AV 

0 

—2umAu  j 
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with  0  =  diag( 0, 0, 0, 0),  05  =  wjvrs,  <76  =  unTq,ct7  =  u>mt7  and  <jg  =  ojmt $.  It  is  sufficient  for  the  proof  if  B 
can  be  shown  to  be  negative  semi-definite.  This  first  leads  to: 

(3.24)  cj,vT6  <  0,  c oNr$  =  umTs,  1  +  unt5  -  umt7  -  0. 


Note  that  we  use  here  the  fact  that  Av  is  symmetrizable  (see  Appendix  B).  Taking  into  account  (3.24),  B 
becomes 


(3.25) 


B  = 


(  2(76 

-<77 

V  -1+^7 


-CTj 

—2un 

0 


—1  +  <77 
0 

—2ujm 


To  ensure  negative  semi-definiteness,  det(B)  <  0  and  therefore 


(3.26) 

Thus 


1  1  \  2  1  „  1 

- 1 - )  g7  “  2 - <77  4-  4(76  4 - 

Wm  OJn  )  OJM  Wm 


<0. 


(3.27)  <  <77  <  a* 

where 

±  _  UN  / (unUm)(  1  -  4(7e {UM  +  Un)) 

7  UM  +  0 ON  y  (un  +  Um)2 

Here  we  note  that  the  condition  that  <76  <  must  be  also  satisfied  in  order  for  <77  to  have  real 

root.  This  yields  the  conditions  in  the  Theorem.  0  Note  that  these  conditions  are  given  independently  of  the 
local  flow  properties.  And  moreover,  the  penalty  parameters  of  each  domain  are  constrained  by  its  adjacent 
domain. 

Remark  2.  Forn  to  be  outward  normal  vector  the  condition  (3.17)  is  now  given  by 


UNTQ  <  0,  UNTQ  +  CJjV/Tg  =0,  1  +  ujnTo  +  umt7  =  0, 


(3.28)  ( —  +  — \  oj2mt 7  +  2r7  +  4 Wjvre  +  —  <  0 

\uM  uN )  uM 

with  the  conditions  (3.13) 

3.2.4.  The  averaging  method.  We  show  in  this  Section  that  the  averaging  method  can  also  be  written 
as  a  penalty  method  with  a  particular  choice  of  the  parameters. 

Euler  Equations  We  start  first  with  the  Euler  equations:  consider  the  following  penalty  method: 

(3.29)  T2,iQ(x,y)[f'~  (q)  -  /'“(?-)], 
where 

(3.30)  f'±={A±qx,A±qy)-ri\go, 

Note  that  the  penalty  terms  use  the  derivative  of  the  fluxes. 
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Theorem  3.3.  If  n  =  73  =  7*2  =  74  —  |  then  the  above  penalty  method  (3.29)  is  equivalent  to  the 

averaging  method  and  is  stable. 

Proof.  We  prove  the  theorem  at  the  interface  x  =  0  with  the  rectangular  domain  and  assume  that 
N  =  M.  If  Ti  =  73  =  and  72  =  n  =  |  then  the  method  becomes 


(3.31) 


<9g7,  _  1  (  dFu\  dG 

dt  Il=°  “  |x=0  _  2  V  dx  +  5a;  )  dy 


and  this  is  obviously  equivalent  to  the  Averaging  Method.  Here  note  that  Following  the 

same  procedure  in  Theorem  3.2,  the  energy  equation  becomes 


2 =  -  (qlAqI  ~  <1Ua<i")  lo  +  ( nq 1  -  r3g77)A+(g7  -  ql1) lo 

(3.32)  +(r2qI  -  r4g77)A~(g7  -  g77) |0. 


Since  n  =  r3  =  r2  =  T4  =  and  g7(0,y,t)  =  g/J(0,  ?/,£),  the  RHS  of  the  above  equation  vanishes  and 
the  energy  is  bounded  by  the  initial  energy.  □ 

The  Navier-Stokes  Equations  :.  The  averaging  method  for  the  N-S  equations  can  be  presented  as 


8 q  8 F  OG  =  OF^  d_G^ 

dt  dx  +  dy  dx  +  dy  + 

ri,3Q(x,y)[fl+(q)-f'+(q-)]  + 
r2,4 Q (x, y) [/'“(g)  -  /'“(«-)]  + 
TsjQix,  y)[Av  ■  a2q  -  A„  •  a2q_]  + 

(3.33)  T6sQ(x,y)[Au  •  dq  —  A„  •  aq_], 


where  a2q  is  the  second  derivative  of  q  in  either  x  or  y  direction. 

Theorem  3.4.  I/tj  =  t3  =  §,  r2  =  t4  =  §,  r5  =  r7  =  §,  and  r6  =  -r8  =  -5^7,  then  the  approximation 
is  continuous  at  the  interface  and  the  scheme  (3.33)  is  stable. 

Proof.  If  n  =  r3  =  |,  r2  =  T4  =  |,  T5  =  7Y  =  5,  and  r6  =  -r8  =  -777,  then  (3.33)  becomes, 

a?7,  _ag77,  _  1  fdF1  dFn\  dG 
dt |l=0  “  dt  u=0  ~  2  V  dx  +  dx  )  dy 

1  (dFl  dFl*\  dGu 

2  \  dx  dx  J  dy 

(3.34)  + — Aj,  *  (dq11  -  dq 7), 

and  this  ensures  the  continuity  of  the  approximation  at  the  interface.  If  the  approximation  is  smooth  enough 
such  that  the  derivative  of  q  is  continuous  at  the  interface  then  this  becomes  the  averaging  method. 

Thus  we  get  for  the  energy: 


— t“E(i)  =  ^  WM1  ~  qUAqn  -  2qIAvqIx  +  2qIIAl/qIx1)  |*=0  - 

fO  f2 

/  qlA^dx  -  /  qIxIAvqIJdx  + 

J —2  JO 

(3.35)  [(rig7  -  r3g77)A+  +  (r2g7  -  r4g77)A” ]  (g7  -  g77)|*=o  + 

[(r5g7  -  T7qn)Al,(qIxx  -  g77)  +  (r6g 7  -  r8g77)A„(g7  -  g77)]  U=o- 


ll 


Since  qI(0,y,t)  =  g/J(0,y,f),  we  have 

2 ~jtE^)  ±  9;([(n  -  t3)A+  +  (r2  -  n)A~U  -  «")  + 

(v>  -  rr)Av(qlxx  -  qii) 

(3.36)  (r6-r8  +  —  )AV(^  -  ^.I))\a=0. 

UN 

Thus  if  ri  =  r3  =  r2  -  r4  =  7-5  —  77-  and  r6  =  -r8  =  the  RIIS  vanishes.  0 

3.2.5.  Adaptive  averaging.  To  ensure  the  stability  of  the  scheme  at  some  particular  collocation  points 
where  the  solution  become  singular  and  unstable,  we  use  the  averaging  method  adaptively  at  selective  grid 
points.  In  particular,  we  switched  from  the  penalty  method  to  the  averaging  when  the  following  criteria  was 
satisfied: 

(3.37) 


or 


(3.38) 


\P-P-\ 

|P  +  P_| 


>Ca 


where  Cave  is  a  non-negative  constant.  Note  that  Cave  =  0  leads  to  the  averaging  method,  whereas  a  large 
Cave  results  in  the  penalty  method.  For  the  value  of  Cave  used  in  this  paper,  we  found  out  that  there  were 
very  few  points  in  which  one  needs  to  switch  from  the  penalty  to  the  averaging  procedure.  Moreover  this 
happened  only  at  very  few  time  steps. 


4.  The  Cavity  System  And  Numerical  Configurations.  In  this  section  we  describe  the  set  up  of 
the  simulations  of  the  recessed  cavity  flameholders  by  the  spectral  multi-domain  technique  presented  above. 
The  main  goal  of  this  experiment  is  to  investigate  how  the  geometry  of  the  aft  wall  affects  the  flame  stability. 

4.1.  Physical  setup.  In  the  SCRAMJet  community,  a  cavity  with  the  length-to-depth  ratio  L/D  < 
7  ~  10  is  usually  categorized  as  an  ’open’  cavity  since  the  upper  shear  layer  re-attaches  at  the  back  face  [2], 
In  this  work,  we  choose  the  L/D  of  the  baseline  cavity  to  be  4  and  thus  the  open  cavity  system  is  considered. 
The  coordinates  of  the  cavity  are  (7cm,  —lcm)  for  the  upper  left  and  (11cm, —2cm)  for  the  right  bottom 
corners  of  cavity.  With  the  length  of  the  neck  of  the  cavity  fixed  to  be  4cm,  we  consider  three  different 
angles  of  the  right  corner  of  the  floor  of  the  cavity  (  60,45  and  30),  we  then  compare  each  one  with  the 
case  of  the  rectangular  aft  wall.  The  fluid  conditions  are  given  as  followings;  the  free  stream  Mach  number 
M  =  1.91,  total  pressure  P  =  2.82 (atm),  total  temperature  T  —  830.6(A)  and  normalized  Reynolds  number 
Re  =  3.9  x  107(l/m).  Note  that  the  Reynolds  number  is  here  normalized  and  has  a  unit  of  1  /[length], 
also  the  Reynolds  number  based  on  the  cavity  dimensions  is  0(1O5).  The  boundary  layer  thickness  scale  is 
6  =  5  x  10“4(m),  and  finally,  the  wall  temperature  is  Tu,  =  460.7835(A).  The  initial  configuration  for  the 
baseline  cavity  system  is  shown  in  Figure  4.1. 

4.2.  Numerical  setup.  We  have  conducted  two  different  experiments  for  each  of  the  following  cases 
(1)  non-reacting  cold  flow  and  (2)  reacting  flow  .  We  use  9  and  17  subdomains  for  both  cases  1  and  2.  For 
the  outflow  conditions  at  the  exit  of  the  system  and  at  the  upper  boundary,  we  mainly  use  a  semi-infinite 
mapping  in  order  to  reduce  the  possible  reflections  at  the  boundaries.  The  characteristic  boundary  conditions 
are  also  applied  and  will  be  discussed  in  the  next  section  and  compared  to  the  mapping.  For  the  case  of  the 
reactive  flows,  the  cavity  was  initially  filled  with  Hydrogen  fuel  with  fuel-to-total  gas  ratio  of  0.5.  The  order 
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Fig.  4.1.  The  initial  configuration  for  the  baseline  cavity  system. 


of  the  polynomial  of  approximation  in  y  direction  in  the  domain  beside  the  wall  is  taken  large  enough  to 
resolve  the  boundary  layer  well.  Finally  the  adaptive  filtering  is  turned  on  if  the  mass  fraction  of  Hydrogen 
and  Oxygen  exceed  the  range  of  -0.09  <  <  1.09,  -0.02  <  fo2  <  0.25  and  the  temperature  exceeds  the 

range  of  300 (A')  <  T  <  3500(A').  As  the  shear  layer  and  the  complex  features  of  the  flows  develop,  the 
adaptivity  criteria  for  applying  the  local  smoothing  is  satisfied  at  some  points.  In  the  calculations,  wre  use 
the  3rd  and  2nd  order  local  filtering  for  the  non-reactive  and  reactive  flows  respectively.  It  turns  out  that 
the  local  smoothing  was  applied  in  very  few  points  at  the  upper  corner  of  the  cavity  wall. 

For  the  adaptive  averaging,  we  use  the  criteria  constant  CaVe  such  that  the  difference  of  the  state  vectors 
(or  pressure)  between  the  two  adjacent  domains  is  less  than  10%.  In  figure  4.2  the  Penalty  Navier-Stokes 
equations  were  considered  for  the  non-reactive  cold  flows.  As  evident  from  the  contours  of  the  density,  the 
approximations  were  well  matched  at  the  interfaces.  Here  the  outer  boundary  was  approximated  by  using  the 
characteristic  conditions  of  the  inviscid  fluxes.  The  adaptive  averaging,  with  the  given  adaptivity  conditions 
above,  took  place  at  only  a  few  points.  The  characteristic  boundary  conditions  using  the  inviscid  fluxes  yield 
good  results  for  both  the  problems  of  the  density  peak  propagation  and  the  non-reactive  cold  flows.  As  in 
figure  1,  we  observe  that  there  exist  penalty  parameters  satisfying  the  stability  conditions  that  may  induce 
reflecting  modes  at  the  interfaces. 


Fig.  4.2.  The  non-reactive  cold  flows  with  the  penalty  Navier-Stokes  equations:  the  density  contours  are  given  in  this 
figure  at  t  =  0.25ms.  17  domains  are  used  and  the  boundaries  of  each  domain  are  shown. 


5.  Results  And  Discussion. 

5.1.  Pressure  history.  Figure  5.1  shows  the  pressure  history  of  the  non-reactive  cold  flows  for  the 
various  angles  of  the  aft  wall  at  two  different  locations  inside  the  cavity,  i.e.  at  the  center,  (x,y)  — 
(8.5 cm,  —1.5 cm),  and  at  the  middle  of  the  floor  (x,y)  =  (8.5 cm,  — 1.9cm). 

These  figures  show  that  the  pressure  fluctuations  in  cavities  with  lower  angle  of  the  aft  are  weaker  than 
in  cavities  with  higher  angles.  It  is  also  shown  that  the  attenuation  of  the  pressure  fluctuations  are  obtained 
both  at  the  center  and  the  middle  of  the  floor  of  the  cavity.  It  is  interesting  to  observe  that  the  patterns  of 
the  pressure  fluctuations  for  a  given  angle  at  different  locations  are  different  depending  on  the  angle.  In  the 
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Fig.  5.1.  Pressure  history  for  non-reactive  flows:  the  left  panel  represents  the  pressure  history  at  the  center  of  the  cavity 
and  the  right  panel  at  the  middle  of  the  floor  of  the  cavity .  Each  panel  shows  the  case  of  90,  60,  46  and  30  degree  cavity  walls 
from  top  to  bottom . 
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Fig.  5.2.  Pressure  history  of  the  non-reactive  flows  with  the  use  of  the  4th  order  filter :  the  left  panel  represents  the 
pressure  history  at  the  center  of  cavity  and  the  right  panel  shows  the  left  panel  in  a  smaller  scale .  Each  panel  shows  the  case 
of  90,  and  30  degree  cavity  walls  from  top  to  bottom.  Note  that  the  scale  of  the  right  panel  is  different  from  the  left. 


case  of  the  30  degree  aft  wall,  the  pressure  fluctuations  are  almost  the  same  at  the  two  locations  considered 
whereas  the  case  of  45  degree  shows  a  difference  in  the  patterns  of  the  pressure  fluctuations  between  the  two 
locations.  The  pressure  fluctuations  at  the  bottom  grows  greater  than  that  at  the  center  after  some  time. 


Fig.  5.3.  Streamlines :  the  left  figure  shows  the  streamlines  at  t  —  1.685ms  for  the  global  filtering  order  7  =  4  and  the 
right  at  t  =  2.38 ms  for  7  =  6. 

Figure  5.2  shows  the  pressure  history  when  the  heavy  global  filter  is  applied  (in  this  case,  the  4th  order 
filter  was  used).  Unlike  the  previous  case  illustrated  in  Figure  5.1,  where  the  6th  order  global  filter  is  used, 
the  pressure  fluctuations  eventually  decay  out  and  a  large  recirculation  zone  is  formed  inside  the  cavity 
without  any  severe  pressure  fluctuations.  Note  that  the  scale  in  the  left  panel  shown  is  the  same  as  in  Figure 
5.1  while  the  right  panel  is  shown  in  a  smaller  scale  for  a  closer  look.  This  figure  shows  that  the  large 
recirculation  zone(s)  formed  inside  the  cavity  obtained  by  the  lower  order  numerical  scheme  is  induced  not 
physically  but  rather  artificially  due  to  the  heavy  numerical  dissipations.  This  is  clearly  shown  in  Figure 
5.3.  In  this  figure  a  large  recirculation  zone  is  observed  -  this  zone  is  formed  earlier  than  this  streamlines 
are  captured  -  when  the  4th  order  filter  is  used(left  figure)  and  an  almost  steady  state  is  already  reached  as 
the  pressure  history  indicates  in  Figure  5.2.  We  find  from  the  numerical  results  that  the  large  recirculation 
is  very  stable  once  it  forms.  This  large  recirculation  and  the  steady  state  solutions  are  not  observed  in  the 
case  of  7  =  6(right).  For  the  case  of  7  =  6  instead  of  the  large  single  recirculation  zone,  smaller  scale  vortex 
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circulations  are  formed  and  they  are  interacting  with  each  other,  never  reaching  the  steady  state  with  time. 
This  result  shows  that  for  these  sensitive  problems,  high  order  accuracy  should  be  used  in  order  to  minimize 
the  effect  of  the  numerical  dissipation. 

Figure  5.4  shows  the  case  of  the  reactive  flows  for  the  90  and  30  degree  aft  walls.  Similar  features  of 
the  pressure  fluctuations  are  shown  as  in  the  non-reactive  flows.  However  the  pressure  fluctuations  are  much 
more  attenuated  for  both  the  90  and  30  degree  walls  than  in  the  non-reactive  cold  flows.  In  the  reactive 
cases  Hydrogen  fuel,  which  was  initially  supplied  inside  the  cavity  was  consumed.  As  time  elapses,  the  fuel 
is  consumed  out  with  the  production  of  the  water  for  these  cases. 


Fig.  5.4.  Pressure  history  for  reactive  flows :  the  left  panel  represents  the  pressure  history  at  the  center  of  cavity  and  the 
right  panel  at  the  middle  of  the  floor  of  cavity.  Each  panel  shows  the  case  of  90  and  30  degree  cavity  walls  from  top  to  bottom. 


Fig.  5.5.  The  water  contour  of  the  reactive  flows:  the  left  the  water  density  contour  is  given  in  the  left  figure  and  its 
streamlines  in  the  right  figure  at  t  —  0.135m$. 

These  results  demonstrate  that  simulations  of  cold  flows  do  not  necessarily  shed  light  on  the  behavior 
of  reactive  flows. 

5.2.  Flow  Helds. 

Non-reactive  cold  flow 

Figure  5.6  shows  the  density  contours  and  streamlines  for  the  90,  60,  45  and  30  degree  walls  at  the 
instant  time  t  =  2.4ms.  As  shown  in  the  figure,  the  shear  layer  is  becoming  weaker  as  the  degree  of  angle 
of  the  aft  wall  and  the  flow  fields  are  becoming  more  regularized  for  the  case  of  the  lower  angle.  And  note 
that  the  density  compression  at  the  corner  of  the  aft  wall  is  also  becoming  weaker  for  the  more  slanted  wall 
cases. 

Figure  5.8  shows  the  streamlines  corresponding  to  the  each  case  of  Figure  5.7.  Note  that  compared  to 
the  non-reactive  cases,  the  shear  layers  are  less  developed  for  the  reactive  cases.  As  the  figures  of  the  pressure 
fluctuation  history  and  Figure  5.8  indicate,  the  shear  layers  are  weak  for  both  the  90  and  the  30  degree  walls 
in  the  reactive  cases. 

Reactive  flow 

Figure  5.7  shows  the  water  contour  inside  the  cavity  for  the  different  angles  at  different  time.  Here  we 
define  the  region  where  the  flames  are  generated  to  be  same  as  the  region  where  the  water  is  produced.  As 
the  Hydrogen  fuel  is  consumed,  the  water  is  produced  and  starts  to  be  expelled  from  the  cavity  to  the  main 
channel.  The  Same-holding  efficiency  is  enhanced  if  the  chemical  radicals  (water  in  this  case)  are  stably 
circulating  and  long  lasting  before  they  are  expelled  from  the  cavity.  Figure  5.7  shows  that  the  lower  angled 
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Fig.  5.6.  The  density  contour  and  the  streamline  of  the  non-reactive  flows:  the  left  column  shows  the  density  contour  for 
90,  60,  4$  and  30  degree  walls  from  top  to  bottom  and  the  right  column  shows  the  corresponding  streamlines  at  t  =  2.43ms. 
The  maximum  contour  level  is  1.8  and  the  minimum  0.5  with  the  level  step  size  50. 

aft  wall  (30  degree  in  this  case)  maintains  more  water  than  the  90  degree  wall  at  a  given  time.  The  figure 
also  shows  that  the  lower  angled  aft  wall  holds  the  flame  (water  in  this  case)  longer  than  the  90  degree  wall 
-  in  the  last  figure  in  Figure  5.7  at  t  =  2.26ms,  the  most  water  is  expelled  and  the  only  the  small  amount 
is  left  in  the  left  corner  while  the  30  degree  wall  cavity  holds  the  water  still  throughout  the  cavity.  These 
results  imply  that  the  flame-holding  efficiency  can  be  increased  by  lowering  the  angle  of  the  aft  wall  of  the 
cavity. 

Appendix  A.  The  similarity  transform  matrices  and  the  eigenvalues  of  the  inviscid  flux 
with  chemical  species. 

Air  model  without  combustion 

First  consider  the  ideal  gas  composed  of  two  chemically  non-reactive  species  (for  the  ideal  mono-atomic 
gas  A  the  diagonal  matrix  and  5,  the  diagonalizer,  were  given  in  [15]).  A  is  given  by 

A  =  diag(U  •  n  +  c,  U  •  ft,  U  -ft,  0  -ft  —  c,  U  -ft,  U  -  ft), 

where  U  =  (u,  v),  ft  =  (nx,ny)  is  an  unit  outward  normal  vector  at  the  interface  and  c  is  a  local  sound  speed. 
For  simplicity  we  assume  that 

p  f  Cp(s)ds  —  P  ~  pCvT , 

Jo 
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t  =  0.175ms 


t  —  0.275 ms 


t  =  0.945ms 


t  =  2.26ms 

Fig.  5.7.  The  water  contour  of  the  reactive  flows:  the  water  density  contours  are  given  in  the  left  figures  for  90  degree  wall 
and  30  degree  wall  in  the  right  figures.  From  top  to  bottom  the  instant  times  t  are  0.175ms,  0.275ms,  0.945ms  and  2.26ms. 
The  maximum  and  minimum  contour  levels  are  0.01  and  0.23  respectively  with  the  number  of  levels  50. 


This  form  is  used  only  in  the  analysis,  as  mentioned  in  Sec.  3.1,  CPi  is  expressed  as  a  4th  order  polynomial 
in  the  temperature  T.  The  nonlinear  expression  of  cpi  makes  it  difficult  to  derive  the  Jacobian  matrices 
of  the  fluxes.  Our  simplifications  is  a  results  of  assuming  small  coefficients  of  the  high  order  terms  of  the 
polynomial.  In  the  actual  simulations  Cv  is  computed  appropriately  using  the  empirical  law  and  assumed 
temperature  independent  at  each  linearization  step.  With  this  assumption  S  is  given  by 
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where  H  =  (E  +  P)fp,Ri  =  R/(MiCv),Rh  =  Rih°2  -  R2h°1,R%,  =  £?=i  =  -!/(*»  +  RD’  the 


t  =  0.175ms 


t  =  0.275ms 


t  —  0.945ms 


t  =  2.26  ms 

Fig.  5.8.  T&e  streamlines  for  the  reactive  flows:  the  streamlines  for  90  degree  wall  are  shown  in  the  left  figures  and  the 
30  degree  wall  in  the  right.  From  top  to  bottom  the  times  t  are  0.175 ms,  0.275 ms,  0.945ms  and  2.51ms. 


tangential  vector  k  =  {-nyynx)  and 


fij  nj 

aij  =  Rv  (Rih^  -  Rjh^y  aji  =  Ru  (Rjh?  -  Rih°j)  • 

Note  that  =  7  -  1  for  the  mono-atomic  ideal  gas  with  7,  the  ratio  between  the  heat  capacities  Cp  and 

Cv. 


Air  model  with  combustion 

Consider  now  the  equations  the  Euler  equations  with  four  reactive  species.  In  this  case  A  and  S  are 
given  by 


A  =  diag(U  •  n  +  c,  U  -  ft,  U  •  n,  U  *  n  -  c,  U  •  n,  U  •  n,  U  •  n,  U  *  n), 
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where  all  the  variables  are  same  as  in  the  two  species  case  except  that 

aijkl  —  tijkl  (ftj  ^2  “b  )#v/#ft, 


Rijkl  —  ^ijkl{Rj  #fc  *b  Hi)cj Rh, 

with 

4 

#/i  “  ^  ^  €ijkiRi{hj  ~  hk  "b  )?  ^5  J?  A;, /  =  1, 2, 3, 4,  j  <  &  <  /, 

i=l 


Cjjki  is  the  permutation  symbol  and  Rv  =  ^=1  .  A  and  5  are  based  on  the  time  dependent  local  spatial 

quantities  at  a  given  time.  is  calculated  at  the  interface  points  at  each  time. 


Appendix  B.  The  symmetrizability  of  the  coefficient  matrices  of  the  Navier-Stokes  equa¬ 
tions  with  chemical  species. 

In  [1]  it  had  been  proven  that  the  coefficient  matrices  of  the  Navier-Stokes  equations  (expressed  in 
the  primitive  form),  of  the  ideal  gas  can  be  simultaneously  symmetrized.  In  [12,  15]  the  same  result  was 
demonstrated  for  the  conservative  form  of  the  equations.  Here  we  show  that  it  is  also  true  for  the  Navier- 
Stokes  equations  of  the  combustible  gas  with  multiple  chemical  species  in  two  dimension. 

Rewrite  the  linearized  Navier-Stokes  equations  (2.2)  in  conservative  form  without  the  chemical  source 
term  as 
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where  A  =  ^,B  =  ^,C  =  |^,D  =  ^  and  E  =  It  is  sufficient  to  consider  the  chemically 

interacting  two  chemical  species.  The  coefficient  matrices  are  given  by 
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where  U2  =  $  ■  U,H  = 

Rvh^,8i  =  -hifi-£-,6i  = 


,0-1 
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=  f  +  E<=1  (fei  -  =  (,-^U2  and  Sji;  =  (TjV?  +  (TfcV2  +  <730  . 

To  find  the  symmetrizer  for  A,B,C,D  and  E,  we  first  consider  the  similarity  transform  matrix  SP  of 
C  such  that 


Sp1CSp  =  Ac, 

where  Ac  is  a  diagonal  matrix  composed  of  the  eigenvalues  of  C.  The  subscript  P  denotes  that  this  matrix 
is  adopted  from  the  parabolic  portion  of  the  equations  [1].  The  diagonal  matrix  Ac  of  C  is  given  by 
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The  diagonalizer  SP  is  composed  of  the  eigenvectors  of  C,  Sp  and  its  inverse  Sp1  are  given  by 
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The  similarity  transform  induced  by  SP,  transforms  the  coefficient  matrices  A,  B,  C,  D  and  E  to 
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where  a  =  Rv0  -  x  and  r)i  =  Rv (Sj  +  hifi)  +  in ■ 

Introducing  a  symmetrizing  diagonal  matrix,  QTQ  such  as 


QrQ  = 


a  0  0  0  0  0 

0  1  0  0  0  0 

0  0  1  0  0  0 

0  0  0  ^0  0 

0  0  0  0  0 


0  0  0  0  0 


21 


we  have  symmetrized  all  the  coefficient  matrices,  i.e. 


QTQSp1ASP  =  (QTQSp1ASp)T, 

QrQSp1BSP  =  (QTQSp1BSp)T, 

QTQSp1CSp  =  (QrQSp1CSp)r, 

QrQSp1DSp  =  (QTQSp1DSp)T, 

QTQSp1ESP  =  Sp1ESp  . 

Appendix  C.  Constants  for  Chemical  Models. 

Here  we  provide  constants  used  in  the  chemical  model  for  the  current  numerical  experiment.  Table  I  gives 
the  constants  used  to  get  the  approximation  of  the  specific  heat  CPi  of  ith  species  in  the  4th  order  polynomial 
of  T,  i.e. 


Cpi  —  (ci  +  T(c2  +  T(c3  +  T(c4  +  c^T))))R/Mi , 
where  jR  is  a  gas  constant,  and  Mi  is  a  molecular  weight  of  ith  species  [16]. 


Table  I 

Coefficients  for  the  approximation  of  the  specific  heat  CPi 


o2 

h2 

H20 

n2 

Ci(l/moJe) 

3.0809 

3.4990 

3.4990 

3.1459 

C2(l/mole) 

0.16962E-2 

-0.18651E-3 

0.14878E-2 

0.99154E-3 

cs(l/mole) 

-0.76334E-6 

0.46064E-6 

0.87544E-7 

-0.22912E-6 

04(1 /mole) 

0.17140E-9 

-0.13157E-9 

-0.11499E-9 

0.12181E-10 

cs(l/mole) 

-0.14116E-13 

0.11679E-13 

0.13495E-13 

0.11024E-14 

Table  II  gives  the  molecular  weight  and  specific  enthalpy  for  each  chemical  species  and  Table  III  gives 
the  reference  dynamic  viscosity,  temperature  constants  T  and  5  in  Wilke’s  law  [21]. 


Table  II 

Molecular  weights  and  specific  enthalpy 


r 

o2 

Ho 

h2o 

n2 

M(l/mole) 

r 

32.000 

2.016 

18.016 

28.016 

h°  (Joule /kg) 

-272918.21  - 

■4280070.46 

-13973684.55 

-302736.23 

Table  III 

Constants  for  Wilke 5 s  law 

o2 

h2 

H20 

n2 

li0{kg/m/sec) 

0.1919E-4 

0.0841  IE-4 

0.1703E-4 

0.1663E-4 

To(K) 

273.111 

273.111 

416.667 

273.111 

S(K) 

138.889 

96.6667 

861.111 

106.667 
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